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Abstract 

Numerical simulation of multiphase compositional flow in fractured porous media, when all the 



species can transfer between the phases, is a real challenge. Despite the broad applications in 
hydrocarbon reservoir engineering and hydrology, a compositional numerical simulator for three- 
phase flow in fractured media has not appeared in the literature, to the best of our knowledge. In 
this work, we present a three-phase fully compositional simulator for fractured media, based on 
higher-order finite element methods. To achieve computational efficiency, we invoke the cross- 
flow equilibrium (CFE) concept between discrete fractures and a small neighborhood in the ma- 
trix blocks. We adopt the mixed hybrid finite element (MHFE) method to approximate convective 

O Darcy fluxes and the pressure equation. This approach is the most natural choice for flow in 

fractured media. The mass balance equations are discretized by the discontinuous Galerkin (DG) 
method, which is perhaps the most efficient approach to capture physical discontinuities in phase 
properties at the matrix-fracture interfaces and at phase boundaries. In this work, we account for 

£^> gravity and Fickian diffusion. The modeling of capillary effects is discussed in a separate paper. 

We present the mathematical framework, using the implicit-pressure-explicit-composition (IM- 
PEC) scheme, which facilitates rigorous thermodynamic stability analyses and the computation of 
^h phase behavior effects to account for transfer of species between the phases. A deceptively simple 

CFL condition is implemented to improve numerical stability and accuracy. We provide six nu- 
merical examples at both small and larger scales and in two and three dimensions, to demonstrate 
powerful features of the formulation. 
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1. Introduction 



Many problems in hydrocarbon reservoir engineering, as well as hydrology, involve the flow 
of multiple distinct phases in fractured porous media. One important example is gas injection in 



Email addresses: j moortgat Or erinst . org (Joachim Moortgat), abbas . firoozabadi@yale . edu (Abbas 
Firoozabadi) 
Preprint submitted to Journal of Computational Physics April 5, 2013 



oil reservoirs that have previously been water flooded. Another example is when gas is injected in 
an oil reservoir and a third hydrocarbon phase develops with intermediate properties to the gas and 
oil phases. When there is significant species exchange between different phases, there is a need for 
multi-phase compositional simulators. A reliable determination of the number of phases, phase 
amounts and phase compositions requires an equation of state (EOS) based phase stability analy- 
sis and three-phase-split computations. Thermodynamic stability analysis is essential to guarantee 
that a phase-split solution corresponds to the lowest Gibbs free energy. For the three-phase- split, 
in particular, there are generally multiple solutions corresponding to local minima. Fully composi- 
tional EOS -based commercial simulators for three hydrocarbon phases are currently not available. 

Compositional simulators commonly use Henry's law or similar correlations to predict the 
C0 2 solubility in water (e.g. |Chang et~aLl ( |1998| )). This is a poor approximation in three-phase 
flow, where it cannot satisfy thermodynamic equilibrium. The C0 2 composition in the aqueous 
phase has to satisfy equality of fugacities of C0 2 in all three phases, which cannot be guaranteed 
by Henry's law. In the presence of C0 2 -rich gas and oil phases, water is not necessarily saturated 
with C0 2 . In our work, we have developed a rigorous EOS-based three-phase thermodynamics 
algorithm. The aqueous phase is modeled by the cubic-plus-association (CPA) EOS, including 
cross-association between C0 2 and water molecules and self-association between water molecules 
( |Li and Firoozabadi||2009| ). At low temperatures (T < 350 K), evaporation of water and the mutual 
solubility between water and hydrocarbon phases can be neglected. To speed up the phase-split 
computation for water-gas-oil mixtures, we only consider C0 2 solubility in water. The CPA-EOS 
reduces to the |Peng and Robinson| ( |1976| ) (PR) EOS for the hydrocarbon phases, when they do not 
contain water. The three-phase compositional model and details of the phase-split computations 
were presented in |Moortgat et al.| ( [2012| ) for unfractured domains. Another compositional model 
for three hydrocarbon phases in unfractured media was considered in |Okuno et al.| ( |20T0| ). In this 
paper, we extend the modeling of both water-oil-gas systems and the flow of three hydrocarbon 
phases to the considerably more complicated fractured porous media. 

Fractured porous media are challenging because of the large range in spatial scales, perme- 
abilities, fluxes and phase properties. Currently, the most efficient compositional simulators are 
based on implicit-pressure-explicit-composition (IMPEC) schemes. The explicit mass transport 
update is constrained by the Courant-Friedrichs-Lewy (CFL) condition on the maximum time- 
step, which is proportional to the size of grid elements, and inversely proportional to the flux 
( |Courant et al4|1928] ). Fractures generally have small apertures, but may allow large fluxes due to 
the high fracture permeability. When fractures are discretized the same way as matrix elements, 
i.e. single-porosity simulations, the resulting CFL condition is exceedingly small and most prob- 
lems are numerically intractable. The most commonly used alternatives are dual-porosity or dual- 
porosity-dual-permeability models, which use two overlapping domains ( |Warren and Root[[l963| ). 
All the flow is through a sugar-cube configuration of fractures, while the matrix only serves as a 
storage medium. The flux between fractures and matrix blocks is computed by so-called transfer 
functions. The dual-porosity and dual-permeability models are adopted in most fractured media 
studies (for various implementations of varying complexity, see for instance |de Swaan| ( |1982| ); 
Coats] ( [T989| ); |Peng et al.| ( [T990| ); [Arana-Ortiz and Rodriguez| ( [T996| ); |Bennion et al.| ( [T999l ); [Cicek 



( [2003] ); [LueFaLl ( [2008] ); |Ramirez etUi] ( [2008] )). The approach is highly efficient, particularly for 
immiscible and single-phase flow, but suffers from severe limitations for multi-phase composi- 
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tional flow, when there is significant species exchange by Fickian diffusion between the fractures 
and the neighboring region in the matrix, gravitational reinfiltration of oil from fractures to matrix 
blocks, or gravitational and viscous instabilities that may cross fractures. Such complex physical 
processes may not rigorously be incorporated in transfer functions. 

We adopt an alternative approach in which fractures are combined with a small fraction of 
the matrix blocks on either side in larger computational elements. The assumption is that a large 
permeability, but small pressure gradient across the fracture-matrix interface results in a transverse 
flux that instantaneously equilibrates the fracture fluid with the fluid immediately next to it in the 
matrix. We call this the cross-flow equilibrium approximation, and denote the combined fracture- 
matrix elements as CF elements. The fluxes across the edges of the CF elements are worked out 
by integrating the appropriate Darcy fluxes for the matrix and fracture contributions. Variations 
of the discrete fracture approach were studied analytically by |Tan and Firoozabadi] ( |1995a|b[ ), 



applied to immiscible water injection in fractured media by |Karimi-Fard and Firoozabadi| ( |2003] ), 



to three-phase black-oil in |Geiger et al.| ( |2009| ); |Fu et al.| ([2005) and to single- and two-phase 



compositional flow in |Hoteit and Firoozabadi| ( [2005 [ |2006[ ). In those papers, and the examples 



presented in this paper for three-phase flow, it is demonstrated that the CF approach provides nearly 
indistinguishable results from fine mesh simulations, but at orders of magnitude lower CPU cost. 
The CF treatment of fractures allows coarser grids, which translates into large CFL time-steps. The 
efficiency of the pressure update is also improved, because of the lower contrast in (CF) fracture 



and matrix permeabilities. In Hoteit and Firoozabadi (2008), the MHFE+DG and CF approach 



was applied to immiscible two-phase flow in fractured media with capillary pressures. In this 
work, we neglect capillarity for simplicity. The additional complications posed by capillarity in 



heterogeneous and fractured domains are presented in a separate paper (Moortgat and Firoozabadi 



2013). 



The paper is organized as follows. In Section [2] the mathematical model is described, fol- 
lowed in Section [3]by the numerical implementation. We discuss the mixed-hybrid-finite-element 
(MHFE) approximation to fluxes and a pressure equation, and construct CF fracture elements 
by appropriately integrating over fracture and matrix fluxes. The discontinuous Galerkin (DG) 
mass transport update, and thermodynamic equilibrium computations are briefly discussed. The 
numerical implementation includes several improvements over earlier work. We extend the CF 
equilibrium discrete fracture model to fully compositional three-phase flow and adopt a more sta- 
ble CFL condition. In Section [4] we provide six numerical examples. Two examples compare 
the CF model to single-porosity simulations, and the other four examples illustrate features of the 
model for both the flow of three hydrocarbon phases, and gas-oil- water systems in fractured two- 
(2D) and three-dimensional (3D) domains. 

2. Mathematical model 

2.1. Mass transport 

We adopt a fractional flow formulation and write the mass- (or species-) transport equation for 
the total molar density of each species i in the three-phase mixture as 

3 c? '• 
0^ + V • U* = F h i=\,...,n c , (1) 

dt 3 



with (p the porosity and c the overall molar density. For each of n c components /, zi is the overall 
molar composition, F t are sink/source terms representing injection and production wells, and U; is 
the total molar flux, which consists of convective phase fluxes d a and diffusive phase fluxes J,- >a : 

U,- = J] (c a Xi,a$a + <f>S a Ji,a) , I = 1, ■ . . , H c . (2) 

a 

The three phases are labeled by a, and for each phase, & a is the Darcy flux, c a the molar density, 
and Xi 9(X the mole fraction of component i. The phases a can either be oil, gas and water, or three 
hydrocarbon phases. The reduction of the diffusive flux by the porosity is a minor improvement 



over the expression in Hoteit and Firoozabadi (2009) to account for the reduced surface available 



for diffusion in porous media. The phase diffusive fluxes are weighed by the saturations for similar 
reasons. An additional reduction may be included to account for tortuosity. 

2.2. Diffusive fluxes 

The diffusive fluxes J it(X are given by 

n c -\ 
<*i,a Ca / , L*ij,a* Xj,ai l A, . . . , Yl c 1, w^/ 

7=1 
n c -\ 



i=\ 



Multi-component Fickian diffusion is modeled by (n c - 1) X (n c - 1) matrices of temperature, 
pressure and composition dependent phase diffusion coefficients, A 7 >, and takes into account 
non-ideality of the fluids. One can easily demonstrate (Hoteit, 1 201 1| ), that mass balance is violated 
when only a single diffusion coefficient, or a diagonal matrix of diffusion coefficient is used, as is 
generally done in commercial simulators. The off-diagonal components represent dragging effects 
and can cause a particular species, in a multi-component mixture, to diffuse from a region of 
low concentration to higher concentration. These effects have been demonstrated experimentally 
for a three-component system ( [Arnold and Toor[ |1967| ). We calculate the diffusion coefficients 



from a unified model for open space diffusion |Leahy-Dios and Firoo zabadi] ( |2007| ). Dispersive 



contributions to D^ a are neglected because they are higher-order in terms of the convective fluxes, 
which are small in the matrix, and lower-dimensional in the fractures. 

2.3. Convective fluxes 

Each of the volumetric convective phase fluxes is given by the corresponding Darcy relation: 

& a = -A a K(Wp-p a gl (4) 

in terms of gravitational acceleration g, permeability tensor K, mobility A a (S a ), saturation S a , and 
mass density p a . A complication of working with individual Darcy phase-fluxes is that Eq. (|4]) 
cannot be inverted in favor of the pressure when a phase may be absent or immobile. To circum- 
vent this issue, and to reduce the system of equations that has to be solved directly, we adopt the 
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fractional flow formalism and write a Darcy relation for the total flux & t . To simplify the notation, 
we adopt the following definitions for phase mobilities A a = k m /fi a , in terms of relative perme- 
abilities k m and viscosities fi a , effective phase mobilities k a = A a K, total (effective) mobilities 
&t - Yia ^a and k t = Yia k<*> an d fractional flow functions f a = A a /A t , and write: 






(5) 



The total (effective) mobility is positive definite, so Eq. ([5]) can be solved for the pressure. As we 
discuss in detail below, we simultaneously solve for the pressure and for d t . After finding the total 
flux d t , we can reconstruct the phase fluxes, independent of the pressure, from: 









(6a) 
(6b) 



2.4. Pressure equation 



We use Acs's method ( |Acs et al.}|1985[|Watts[|1986| ) to compute the pressure field from: 

dp 



0G-^ + ^v ; <V-U ( -F ; ) = O, 



(7) 



where C t and v,- are, respectively, the total compressibility and total partial molar volumes of the 



three-phase mixture. Expressions for both variables are derived in Appendix C in (Moortgatet al. 



2012). Rock compressibility may also be included in C t . 



2.5. Phase compositions and molar fractions 

Phase compositions x^ a and phase molar fractions m a are derived from the non-linear set 
of equations that guarantee equality of fugacities of each component i in all three phases (a = 
0fi,0f2,ar3), as required by local thermodynamic equilibrium. Computationally, the natural loga- 
rithm of the equilibrium ratios K ia is more robust (H augen et al.[|20lT] ). Selecting one reference 
phase, say a^, we have two sets of equilibrium ratios K^ ai = x^ ai /x^ a3 and K^ a2 = x^ a2 /x^ a2 
satisfying the equilibrium criteria: 



hi K 



i,a\ 



= ln<^ >3 -ln^, and 



(8a) 
(8b) 



ln^ >2 = ln<^ >3 -ln<^ >2 , 
in terms of the fugacity coefficients <p ita . Eq. ([8]) is supplemented by the constraint relations: 

Zt = 2_j m <* x i,<*> i = 1, • • • , n c - 1, (9a) 

a 
i i i i 
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An equation of state has to be specified to describe the three phases and derived quantities, such 
as saturations, molar and mass densities, compressibilities and partial molar volumes. We model 
the aqueous phase with a cubic-plus-association EOS that takes into account cross-association 
between water and C0 2 molecules, and self-association of water. In the absence of water, the 
CPA-EOS reduces to the Peng-Robinson EOS, which we use for pure hydrocarbon phases. We 



refer the reader to (Li and Firoozabadi, 2009; Moortgat et al. , 2012) for details 



2.6. Boundary and initial conditions 

To complete the description of the physical model, we prescribe initial and boundary condi- 
tions. The initial condition consists of the overall composition and pressure field throughout the 
domain. The boundaries are described by non-overlapping Dirichlet and Neumann conditions: 
we consider impermeable boundaries, except in production wells where we have either a con- 
stant pressure or production rate. Injection wells are placed inside the domain as source terms. 
Production wells can also be described as sink terms, with impermeable boundaries everywhere. 

3. Numerical model 

3.1. Mixed hybrid finite element method 
3.1.1. Expansion of convective Darcy flux 

In the MHFE method, the convective and gravitational fluxes are decomposed into their normal 
components across the edges E of each computational matrix element K as: 

d(t,x) = ^j fe^(0w^(x), and g(x) = J^ q K K g E ^ KE (x), (10) 

EedK EedK 

where x = (x, y), dK is the boundary of element K, and w K ^ E are the lowest-order Raviart- Thomas 
basis vector fields (Raviart an d Thomas! |1977| ). These vector functions satisfy the properties 



v?k,e ' ^k,e' = — , and V • vr KtE = — , (11) 

Ae Vk 

where A E is the length/area of edge/face E, V K is the area/volume of element K. The MHFE 
weak form of Eq. ([5]) is obtained by multiplying by vr KtE and integrating over each element 
K. The pressure gradient term is partially integrated and Gauss' theorem is used, such that we 
have one volume integral over the pressure, and one surface integral. We define p K = f p and 

L K P - Tje f E P = He tpK,E, which are the averaged pressure in a matrix element, and the aver- 
aged pressures along the element edges/faces. We refer to the latter as pressure traces. The MHFE 
approximation to Darcy 's law can then be written as: 



<lK,E 



= 0k,ePk ~ ^ PK,E,E>tp K ,E> ~ 7k,e, E £ dK. (12) 



E'edK 



The coefficients KiE , J3 k ,e,e' an ^ 7k,e are defined in Appendix |A| 

In two dimensions, fracture elements are initially treated as ID computational elements, such 
that volume integrals reduce to line integrals over fracture elements /, multiplied by the fracture 
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aperture e. In three dimensions, the fractures are 2D planes with an e width in the third dimension. 
Similar to the definitions above, we denote the average pressure in a fracture element by p f , the 
pressures at the end-points or end-edges e of df by tp f , and the size of fracture element / by \f\. 
The MHFE expression for fracture convective fluxes is: 

q/,e = Of,ePf ~ J] Pf^'tPW ~ 7fa> e G d f' ( 13 ) 

e'edf 

However, as discussed in the next section, there is no need to evaluate Eq. ( [T3] ) explicitly. 

3.1.2. Cross-flow equilibrium approximation 

An explicit treatment of fractures, referred to as single-porosity models, is not computationally 
feasible. The CFL constraint on the time-step scales with At oc min(V^/^, Vf/qf). For fracture 
elements, the fracture flux q f may be high, while the volume Vf of a fracture element is generally 
exceedingly small, particularly for fracture intersections. To overcome this limitation, we note 
that the pressure field is continuous, so the pressure in a fracture element is close to the pressure a 
small distance away in the matrix. We now represent a fracture element, together with two small 
matrix elements on both sides, as one computational element. This is illustrated in Figure [T] for a 
rectangular 2D mesh. For this larger, combined element, we assign one averaged pressure p K and 
4 (in 2D) or 6 (in 3D) pressure traces tp K , E . 

Fluxes through edges that are intersected by a fracture (top and bottom edges in Figure [T]) are 
computed by properly integrating over both the matrix and fracture fluxes inside the combined 
element. We now use K for an element that may contain a fracture, k for the matrix portion of 
element K and E for an edge that may be intersected by a fracture and write for the total fracture 
plus matrix flux q KyE = q kiE + q fe , with 

(lK,E = (Ok,E + Of,e)PK ~ /_j {P k > E > E ' +Pf,e,e'<5e,E<5e>,E>) tp k ,E> ~ 7k,E ~ 7/,A,E, (14) 

E'edK 

where 5 e , E is 1 when an edge E is intersected by a fracture element end e and zero otherwise. 
Eq. ( [T4| ) is similar to the two-phase expressions in |Hoteit and Firoozabadi| ( |2006| ). However, in 



Appendix [A] we demonstrate that Eq. ( [T4| ) can be more elegantly reduced to Eq. ( [T2| ), with the 



coefficients k ,e 9 @k,e,e' and y KiE evaluated in terms of a weighted total effective mobility k^ ff 
across edge E: 

kf = k tJ s/A E + k, m (l - e/A E ), (15) 

where s is the area of the fracture intersection with E, and the subscripts m and / denote matrix 
and fracture properties, respectively (i.e. \ m is the total effective mobility in the matrix, as defined 
above Eq. ([5])). We emphasize that we allow different relative permeabilities in the fracture and 
matrix portions of CF elements. 

Fluxes in the transverse direction (left and right edges in the figure) are matrix fluxes, computed 
from Eq. ( [T2| ). The fracture-matrix flux inside the element is accounted for by the assumption that 
at the fracture-matrix interface there is a large permeability, but a small pressure gradient, which 
results in a Darcy fracture-matrix flux. The assumption is that this flux instantaneously equili- 
brates, or mixes, the fluid in the fracture with the fluid in the small neighboring matrix elements 
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(small with respect to the full matrix block). When the MHFE method is combined with a FD mass 
transport update, this means that the average czt for the combined element is updated. For the com- 
bination of MHFE with a higher-order DG method, czi at the nodes or edges is updated. The mass 
transport update, Eq. ([T]), is unaltered from its implementation for homogeneous media. We pro- 
vide the phase compositions for the combined fracture-matrix element and the summed fluxes, 
and update the overall molar species densities czu We refer to this approach as the cross-flow (CF) 
equilibrium model and will refer to the combined fracture-matrix elements as CF elements. 

In a single-porosity simulation, the CFL condition could be determined by fracture-intersection 
elements with an area of, say 1 mm 2 , and fracture fluxes that scale with the high fracture perme- 
ability. In the CF approach, when we have, 10 m x 10 m matrix blocks, we can use CF elements 
with a width of several cm or more. This results in a CFL constraint on the time-step that is several 
orders of magnitude larger than for a single-porosity model. A second reason that single-porosity 
models are computationally expensive, is that the system of equations in the pressure update is 
ill-conditioned due to the high permeability contrast between the fracture and matrix elements. 
When we use the averaged CF elements, the pressure update (discussed below) is considerably 
more efficient. 

We emphasize, that the matrix blocks may be discretized by any number of grid-cells, such 
that we can resolve potential gravitational or viscous fingers in the matrix. The discretization of 
the matrix blocks is particularly important to model diffusion. The diffusive fluxes are weighed 
by the phase saturations (Eq. ([2])), and diffusion only occurs within a given phase. When gas is 
injected in fractured porous media, all the hydrocarbons in the fractures may evaporate into the 
gas phase, before a gas phase has developed in the neighboring matrix blocks. Because Fickian 
diffusion only occurs within a phase, one cannot self-consistently compute a diffusive flux from the 
fractures to the matrix blocks. In reality, the gas and oil at the fracture-matrix interface are in local 
thermodynamic equilibrium, and dissolution of gas and evaporation of oil occur through Fickian 
diffusion at the interface. This numerical issue is particularly problematic in single- and dual- 
porosity models. In the CF approach the problem is alleviated, because the gas in the fractures 
is mixed with matrix oil and the CF elements may remain in two-phase, particularly in large- 
scale simulations where breakthrough is avoided. When a sufficient amount of light species has 
accumulated in the neighboring matrix elements, a gas phase may form and diffusion can occur 
from fracture to matrix in both phases at a high rate. In this fashion, the light species may diffuse 
element by element into the matrix blocks, while heavier species diffuse towards the fractures. 

Gravitational effects, such as fingering and re-infiltration, where oil drains from a matrix block 
into a fracture, and then drains from the fracture into another neighboring matrix block, are mod- 
eled without special treatment. The CF model is not restricted to sugar-cube fracture configura- 
tions, but can be applied to any configuration of discrete fractures in structured or unstructured 
grids. These features are difficult to incorporate in the dual-porosity model. 

3.1.3. Pressure equation 

The discretization of the pressure equation, Eq. ([7]), is also greatly simplified by the defini- 
tion of the weighted effective mobility in CF elements. For the convective Darcy flux through a 
fracture-intersected edge/face E, one can easily see that the total phase flux through the matrix and 
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the fracture portions of the CF element reduces to 

q a ,K,E = -kf (Vp - p a q 8 K , E )- (16) 

Similarly, we can use Eq. (J6J) with k a replaced with k^ ff and f a with fjf . 

We expand the diffusive fluxes, similar to the convective fluxes in Eq. ( [T0| ) as: 

(f>SJ Ua (t,x) = Yj *£*e(')w*e(x). (17) 

EedK 

The diffusion term in Eq. (|7J) (with Eq. ([2])) then reduces to 2/ Yje vrtfaKE anc ^ can ^ e combined 
with the source/sink term in F t . For brevity, we move these terms to the right-hand- side of the 
equations, denoted as r.h.s. and define £ = C t ypfVf + 4> m V k ). The integral form of Eq. M is: 

r\ Tl c x-» 

£~5T + yi p * I ( m i,K$t,K ~ Sijd ' n K = r.h.s., (18) 

where in CF elements #>^ is the total flux integrated over both fracture and matrix portions of 
element K, and m UK and s^ are defined in Appendix [Bj In other words, the discretized pressure 
equation has the same form for fracture-containing CF elements as for matrix elements, i.e. is 
the same for fractured and unfractured domains, when written in terms of the weighted effective 
mobilities and £. This formulations is more straightforward than the earlier implementation for 
two-phase flow in poteit and Firoozabadi| ( [2006| ). 

Expanding the fluxes as in Eq. ( [T0| ) and carrying out the integrations using Eq. ( fTT] ), we find 

£-7jT + 2 ^ Jj ( mi > K > E V K > E " S i>K#) = r ' h ' S - ( 19 ) 

i=l EedK 

We eliminate the fluxes by Eq. ( [T2| ) to obtain the spatial discretization of the pressures: 

£-tjT + &kPk ~ ^PKjtpzj ~7k = r.h.s. (20) 



EeK 



For the temporal discretization we use the backward (implicit) Euler method. The fully discretized 
pressure equation, with the r.h.s. re-instated, becomes: 



pr' = ^ A ' l( 



EeK i \ E a )) 



with all the coefficients evaluated at the previous time- step. 

3.1.4. Assembly of global matrix for the pressure update 

Eqs. ( [T2| ) and ( [2T] ) are for individual elements and edges. To construct the global system of 
equations to solve, we assume flux continuity across element edges: 

q K ,E + qK>,E = 0, for E = KHK\ (22) 
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(note that fluxes are defined with respect to the normal to edge E in element K). Collecting the 
terms for each edge in the domain, we obtain the matrix system 

R r P - MT P = I (23) 

For N K grid elements and N E edges, P is the A^-size vector of element averaged pressures, T P is 
a A^-size vector of the pressure traces at all the edges. The matrices R and M, and vector I are 
defined in Appendix [C] 

Similarly, we assume pressure continuity, such that at each edge E 

tp K ,E = tp K >,E, for E = KHK\ (24) 

and construct the matrix system (with matrix definitions given in Appendix [C]): 

DP - RT P = G. (25) 



Eqs. ( [23] ) and ( [25] ) can be combined in one large system that solves simultaneously for the pressures 



and pressure traces. A more efficient approach presents itself by the fact that D is diagonal. By 



multiplying Eq. (25 ) by D *, we can eliminate P from Eq. (23) and obtain a system for T P alone: 

(M - R r D _1 R)Tp = R r D _1 G - I. (26) 

Eq. ( [26] ) is the system that is solved in the pressure update. The matrix that needs to be inverted has 
dimensions N E x N E , but is sparse. On a structured 2D (3D) grid, each (non-boundary) row/edge 
has 7 (11) non-zero elements, because tp K ^ E depends on edge E and the other 3 (5) edges in 
each of the two neighboring elements. After updating T p , P is found through inexpensive back- 
substitution in Eq. ( [25] ). The total flux is found by substituting the updated p K and tp K in Eq. ( [T4| ). 



The system Eq. ( [26] ) is larger than for a FD method, but the considerable advantage is that 
we simultaneously solve for the element pressures (as in FD), the pressures on the edges (which 
is advantageous for heterogeneous and fractured domains), and for a continuous velocity field 
throughout the domain; all with the same order of convergence. 

3.2. Discontinuous Galerkin mass transport update 

All higher-order methods approximate the mass transport update, Eq. ([]]), by multiple degrees 
of freedom for the overall and phase compositions and molar densities. The DG method has the 
additional advantage that the variables can be discontinuous across edges. One of the benefits is 
that different orders of approximation can be used in different elements. In our work, we use a 
bilinear (trilinear) approximation on all 2D (3D) structured elements and a linear approximation on 
unstructured triangular grids ([Moortgat and Firoozabadi[[2010|). Away from strong compositional 



gradients, we can use a lower-order approximation. Our main interest in the DG method is the 
simulation of heterogeneous and fractured porous media. At the boundaries between regions of 
different permeabilities (such as fractures and layers), the phase properties may exhibit strong 
discontinuities. Continuous higher-order methods, such as FD and finite volume, may be used in 
homogeneous domains, but are a less natural choice to approximate the inherently discontinuous 
phase properties in fractured reservoirs. 
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The higher-order mass transport update reduces numerical dispersion and converges to the 
exact solution at a higher rate. Alternatively, this means that a given FD result can be reproduced 
on a significantly coarser grid, and at a correspondingly higher CPU efficiency. Combined with the 
accurate velocity field from the MHFE method, this approach has been shown to result in orders of 
magnitude improvement in CPU times (Ho teit and Firoozabadi[|2005j|2006[|Moortgat et al.[|201 1) 
Moortgat and Firoozabadi[|2010^|Moortgat et al. , 2012 ). A convergence analysis in |Moortgatet al.| 
( 2013 ) demonstrates twice the convergence rate for the DG mass transport update as compared to 
a FD approach. 

As was mentioned above, the implementation of the DG method is identical to that in homo- 
geneous domains, which was presented for three-phase flow in |Moortgat et al.|(2011[ ) for 2D and 
Shah raeeni et al.| ( |2013| ) in 3D, and will not be repeated here. Phase properties are updated at 
either the edge-centers or the nodes from the total (fracture plus matrix) fluxes through each of 
the edges. The accuracy is further improved, because we have the pressures at the edges from the 
MHFE update. To avoid spurious oscillations that may occur in higher-order methods, we use the 
same slope limiter as in the papers cited above. 



3.3. Phase behavior 

We have developed a phase splitting package that can model both three hydrocarbon phases, 
with transfer of all species between the three phases, and systems in which one of the phases is wa- 
ter. In the latter case, we neglect the mutual solubility between water and hydrocarbons and water 
evaporation, and only allow C0 2 solubility in water. These are reasonable assumptions for prob- 
lems where the temperature is below 350 K, and result in considerable computational advantage. 
In the update of phase compositions, first a stability analysis is performed, corresponding to the 
minimum Gibbs free energy. Then, two- or three-phase- split computations are carried out. When 
initial guesses are not available, the phase-split routine first performs a number of successive- 
substitution-iterations (SSI) to obtain a good enough initial guess to switch to the fast-converging 
Newton method. Generally, the Newton method, based on the natural logarithm of the equilibrium 
ratios (Eq. ([8])) only needs one or two iterations to converge. When initial guesses are available 
from the previous time-step, Newton's method is attempted first, without a stability analysis. Var- 
ious optimizations have resulted in a highly efficient algorithm. The details are provided in an 



earlier paper on three-phase flow in homogeneous media (Moortgat et al., 2012). 



3.4. Relative permeability and viscosity 

We adopt Stone I relative permeabilities ( [Stone! [T9701 11973| ). In the numerical examples, we 



will denote the residual saturations for water-oil-gas mixtures by S" w for water, S" ow for oil-to- 
water, S^ a for oil-to-gas, and S® for gas. Similarly, the end-point relative permeabilities are /^ w , 



e 



rog 

/^ og , and k® g , and the powers are n w , n ow , n og , and n g . For mixtures of three hydrocarbon 
phases, we will use the same notation for simplicity, where the w subscript refers to the third 
hydrocarbon phase. 

We use either the LBC ( |Lohrenz et al.[ |1964[ ) or PC ( |Pedersen et al.[ |1984| ) viscosity corre- 
lations. Near the critical point, the LBC correlation may perform poorly and the PC correlation, 
which does not require phase identification, is an improvement. 
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3.5. CFL condition 

Earlier work on higher-order IMPEC modeling of two-phase compositional flow in fractured 



media ( Hoteit and Firoozabadi , 2005 , 2006 , 2009 ) used the CFL condition for the convective fluxes 



(AOconv < min I ^L -) . (27) 

K \Lee8K La Ha,K,E\) 

However, this condition only applies to immiscible flow and does not guarantee that the total 
number of moles of any species cannot flow out of any grid element in one time-step. Moreover, 
the summation over the absolute values of the fluxes may result in an overly restrictive time-step 
when some of the fluxes flow into an element. 

We have implemented a different CFL condition for the compositional fractional flow formu- 
lation, equivalent to the suggestion in |Coats| ( |2003| ), in terms of the convective fluxes: 



(AOconv < min I J^kZ4x _ j ? y ~^ > ^ (2g) 

hK \ LEedK La c a,K,E^i,a,K,E^la,K,E / 

where c a ^ E and x Ua ^ E are, respectively, the phase molar density and composition, evaluated on 
edge E in the higher-order DG mass transport update. c K and zi,k are the element averaged molar 



density of the mixture and overall composition in element K. By using condition Eq. 28 , the 
stability of the algorithm is improved considerably, in particular for fractured media. 

When Fickian diffusion is included, we add the outgoing diffusive fluxes to the denominator 
in Eq. ( [28] ) and check an additional criteria in terms of the maximum eigenvalue of the matrix of 
diffusion coefficients. Denoting the maximum eigenvalue of each of the matrices of phase diffusion 
coefficients by Adiff,^, and A diff ^ = max Qf (A diff?Qf ^), we define 

(AOdiff < minf-^L), and (29) 

K \Adiff,W 

At < min((A0conv, (AOdiff) (30) 

However, most problems of interest are convection dominated. In homogeneous media, the con- 
vective fluxes are generally larger than the diffusive fluxes, even at low injection rates. Fickian 
diffusion is most pronounced in fractured media, where steep compositional gradients may exist 
between fractures and matrix blocks. However, the CFL condition is usually determined by the 
largest convective flux inside small fracture elements. Compared to FD models, the CFL constraint 
is alleviated, because the higher-order DG method allows the use of coarser grids elements. 

The steps in the full simulator algorithm are similar to other IMPEC codes, and are outlined in 



more detail in Moortgat and Firoozabadi (2010); Shahraeeni et al. (2013) 



4. Numerical examples 

We present six numerical examples to illustrate the strengths of our discrete fracture three- 
phase flow model. In Examples 1 and 5, we compare cross-flow equilibrium results to single- 
porosity simulations. In these example, we neglect Fickian diffusion, because single-porosity 

12 



simulations have a numerical issue in computing diffusion in fractured media, as discussed in 
Section [3X21 

In Examples 2 and 3, we consider a typical oil recovery scenario for fractured porous media 
and account for diffusion. The domain is depleted, followed by water flooding and then enhanced 
oil recovery by C0 2 injection. Example 3 considers a high matrix permeability, such that gravi- 
tational fingers may develop throughout the domain. The flow of three hydrocarbon phases in a 
larger fractured domain is illustrated in Example 4. Examples 1 - 5 are for 2D domains. In the 
last example, we consider C0 2 injection into a complex 3D domain with a number of discrete 
horizontal and vertical fractures. 

4.1. Example 1: Comparison of discrete CF model to single porosity simulation 

We test the CF model by comparing to a single-porosity simulation for water flooding, followed 
by CO2 injection in a 2 m x 10m column with four matrix blocks. Convergence of the CF results 
is verified by performing simulations on a coarse 1 1 x 57 element Grid 1, and a finer 21 x 105 Grid 
2. The domain, grids and locations of fractures and wells are indicated in Figure [2} To make the 
single-porosity simulation computationally feasible, we assume relatively wide fractures of 5 mm. 
For the cross-flow simulation we use a 10 times larger width of 5 cm. The fracture permeability is 
4 d, the matrix permeability is 4 md and matrix porosity is 15%. 

The column is initial saturated with a light oil, with composition and EOS -parameters given in 
Table [TJ The density and viscosity of the oil, water and C0 2 at the initial condition of T = 350 K, 
Po = 300 bar are given in Table [2} Both water and C0 2 densities are higher than the oil density, 
and both are injected from the bottom fracture. Production is at constant pressure from the top 
fracture. The relative permeability parameters for all the examples are listed in Table [3} 

In the simulations, we first inject one pore volume (PV) of water and then 250% PV of C0 2 at 
2% PV/day. The purpose is to verify that we can obtain the same results with our CF model as with 
a single-porosity simulation, but at high CPU efficiency. Figure [3] shows the oil saturation at the 
end of water flooding for a single-porosity simulation, and for CF simulations on Grids 1 and 2. 
The oil saturations at the end of the simulations are given in Figure [4j We find excellent agreement 
between the single-porosity and CF results, and observe that the CF results have mostly converged 
on the coarser Grid 1 . The CF simulation is 26 times faster than the single-porosity simulation 
(1 .23 hr versus 32 hrs on Grid 1). The CFL condition is determined by the grid elements containing 
fracture intersections, which are about 24 times larger for the CF discretization (CFL ~ 1 hr) than 
for the single-porosity elements (CFL ~ 2.5 mins). The lower contrast in permeability between 
CF and matrix elements also results in a more efficient pressure update. Note that for a smaller 
fracture aperture the difference in CPU time would be orders of magnitude, and the single-porosity 
simulation would not be feasible, even for this small problem. 

Figure [5] shows the oil recovery for both simulations. The results are close, but the recovery 
from the CF simulation is about 2.2% higher. The reason is that when the gas in the fractures is 
mixed with a small amount of matrix fluid, the light oil evaporates and is quickly recovered in this 
small-scale example, where breakthrough occurs early. When we subtract the small pore volume 
of oil included in the CF elements, the recoveries from CF and single-porosity simulations are in 
perfect agreement. Also shown is that the CF results on Grids 1 and 2 have converged in terms of 
oil recovery predictions (the CPU time for the CF simulation on Grid 2 is 3 hrs). 
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Further verification of the CF model for single- and two-phase flow was made in |Hoteit and 



|Firoozabadi] ( |2005[ |2006| ). Now that we have confidence in the CF model, we proceed to a similar 



problem, but on a large scale that is computationally too expensive for a single-porosity simulation. 

4.2. Example 2: Depletion, water flooding and C0 2 injection in fracture porous media 

We consider a 50 m x 50 m domain with 5 m x 5 m matrix blocks. Again, we verify 
convergence on three different mesh refinements. The grids, fractures and wells are shown in 
Figure |6j and results will be illustrated on Grid 3, unless stated otherwise. The matrix-blocks 
have a porosity of 20% and permeability of 10 md. The fractures have an aperture of 1 mm and 
permeability of 25 d, and the width of the CF fracture elements is 20 cm. The domain is initially 
saturated with oil, with composition and EOS parameters given in Table [4j The temperature is 
350 K and at this temperature, the fluid has a bubble point pressure of 338 bar. 

The simulation is started at an initial pressure of 350 bar at the bottom of the domain. At this 
pressure, the entire domain is in single-phase. In the first recovery phase, we deplete the domain 
at a constant rate of 5% PV/yr (computed at the initial pressure) for one year, until the pressure 
has dropped to 300 bar. Below the bubble point, a gas cap develops in the top of the domain. 
Table [5] provides the densities and viscosities for water, oil, C0 2 , and liberated gas at T = 350 K 
and p = 300 and 350 bar. The water and C0 2 densities are higher than the oil density, so both are 
injected from the bottom. 

Figure [7] shows the gas and oil saturations at the end of the depletion stage. The gas, liberated 
at the reduced pressure, has segregated to the top through the high-permeability fractures. During 
secondary recovery, 50% PV of water is injected at a constant rate of 5% PV/yr from the bottom 
well and production is at a constant pressure from the top. Figure [8] shows the water saturation 
at 5% and 50% PVI. Because of the residual oil saturation to water, breakthrough has occurred 
and further water injection is inefficient. To recover the residual oil, we consider tertiary recovery 
by C0 2 injection for 20 years at 5% PV/yr. The overall and three phase compositions of C0 2 at 
one PV of C0 2 -injection are shown in Figure [9| In particular, we note the C0 2 dissolution in the 
aqueous phase in the large three-phase region. C0 2 has a solubility of 2.5 mol%, or about 6 wt% 
at the given pressure and temperature. The C0 2 that dissolves in water is lost to the oil sweep, but 
results in swelling of the aqueous phase by ~ 5%. At the same time, C0 2 has a high solubility in 
oil which leads to swelling of the oil phase as well. 

Figure [10] shows the oil recovery from depletion, water flooding and C0 2 injection. Oil re- 
covery from depletion is 5%. Water flooding produces another 30% recovery and the enhanced 



oil recovery by C0 2 injection achieves a final recovery to 60%. Figure [TO] also shows that the 
DG results have converged even on the coarsest mesh, in which the matrix blocks are discretized 
by only three elements in each direction. The CPU times for this simulation are 10, 19 and 34 
hrs for Grid 1, 2 and 3, respectively. Because the domain is fractured, phase states and composi- 
tions vary wildly throughout the domain and 95% of the CPU time is consumed by the phase split 
computations. As an illustration of the efficiency of the CF model itself: without the phase- split 
calculations, the CPU time on the finest mesh is only 2.5 hrs. The computational cost of three- 
phase split computations in fractured media motivates continued efforts in improving the efficiency 
of these algorithms. 
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4.3. Example 3: Gravitational fingering in high-permeability porous media 

We repeat the example above but increase the matrix permeability by a factor 10, and consider 
a temperature of 400 K, such that the C0 2 is lighter than the oil (Table [5]). After depletion and 
water flooding, C0 2 is injected from the top well, and production is at constant pressure from the 
bottom. The simulations are carried out on Grid 2 in Figure [6} and account for Fickian diffusion. 

When C0 2 dissolves in water, it increases the density of the aqueous phase by about 1%. 
Because C0 2 is injected on top of the previously injected water, the density increase first occurs in 
the top. Even a density change this small may be gravitationally unstable and trigger a fingering 
instability. The dissolution of C0 2 in oil results in a similar density increase in single-phase, and a 
higher density increase when light components evaporate from the oil into the C0 2 -rich gas phase. 
When the matrix permeability is low, the gravitational fingers do not develop and may be stabilized 
by Fickian diffusion. At higher permeability, the fingering speed scales linearly with permeability. 

Figure [TT] shows the C0 2 composition in the aqueous phase and resulting density increase after 
depletion, water flooding and injection of 30% PVI C0 2 . Pronounced gravitational fingers have 
developed that span multiple matrix blocks and fractures. Complicated flow patterns like this may 
not be studied with dual-porosity models. In FD simulations, gravitational and viscous instabilities 
are often suppressed by numerical dispersion, unless unreasonably fine meshes are used. This 
example illustrates many of the powerful features or our model: the density increase of water, 
predicted by the CPA-EOS in three-phase flow, the high accuracy and low numerical dispersion of 
our higher-order finite element methods, and the efficient modeling of fractures (CPU time of 5.3 
hrs). 

4.4. Example 4: Flow of three hydrocarbon phases in fractured porous media 

We consider the flow of three hydrocarbon phases in a fractured domain, with transfer of all 
species between the three phases, but neglecting diffusion. The domain, fracture network and 
aperture, and porosity are the same as Grid 2 in the previous two examples. The matrix and 
fracture permeabilities are 4 md and 40 d, respectively. The domain is initially saturated with the 



North Ward Estes oil, and the fluid composition and EOS -parameters can be found in Moortgat 



etaE1 ( [20T2| ); |Khan et al.| ( [T992| ); |Okuno et al.| ( [20T0| ). The temperature is 301 K, the initial pressure 



at the bottom is 75 bar, and the relative permeability data are provided in Table [3j 

We inject a mixture of 95 mol% C0 2 and 5 mol% methane from the top at 10% PV/yr for 



10 years, and produce at a constant pressure from the bottom. Figure [T2| shows the overall and 
three phase molar compositions of C0 2 at 100% PVI. In the region between the lightest C0 2 -rich 
gas phase in the top and densest oil phase in the bottom, a large intermediate C0 2 -rich phase has 
developed (denoted by C0 2?z ). This third phase has a low residual saturation (5%) and is readily 
recovered, while the oil phase is stripped from some of its lighter components and remains as 
a denser more viscous residual oil. Such flow properties can not be captured by a two-phase 
compositional simulator. Moreover, we find, in line with remarks in |Okuno et al.| ( |2010| ), that 
modeling an intrinsically three-phase problem with a two-phase simulator may result in erratic 
behavior and crashes. In three-phase regions, the phase compositions obtained from a two-phase- 
split computation do not correspond to the lowest Gibbs free energy (which is the three-phase 
result). The derived phase properties may be unstable to small perturbations. 
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The oil recovery for this example is provided in Figure [T3| The recovery is not intended to 
be representative for a field scale project, in which C0 2 injection should only be considered for a 
shorter duration, and certainly not long beyond breakthrough. 

The total CPU time for this example is 6 hr, the average CFL time- step is 5 hr, the CPU time 
per time-step is one second, with 97% of the computation time spent on the phase-split com- 
putations. As was mentioned above, the cost of phase-split computations in fractured media is 
considerably higher than for homogeneous domains, because phase boundaries occur throughout 
the domain. Away from phase boundaries, most phase-split computations can be avoided by cri- 
teria to determine whether an element was and remains in single-phase. Stability analyses can 
be skipped in two- and three-phase regions when initial guesses are available from the previous 
time-step. In fractured media, compositions vary significantly throughout the domain, and full 
stability and phase-split computations have to be carried out in most elements. For this exam- 
ple, 14% of three-phase split computations were carried out with initial guesses from the previous 
time-step, avoiding the stability analysis and two-phase split, and only using a couple of Newton 
iterations. Another 15% of phase- split calculations were avoided by determining which elements 
were and remained in single-phase. As an indication of the efficiency of the MHFE+DG and CFE 
fracture model, the CPU time for this example, with the phase-split computations subtracted, is 
only 8 mins. The large fraction of CPU time spent on the phase-split computations is of course 
partly due to the small mesh size of only 2401 elements, for which the pressure and mass transport 
updates are extremely fast. We also note that we compute phase compositions to an accuracy of 
10" 10 , and the final mass balance error for each of the components is of the order of 10" 15 . 

4.5. Example 5: Comparison of CF and fine mesh simulations for three hydrocarbon phases 

We compare briefly the CF model to a single-porosity simulation for a three hydrocarbon 
phase fully compositional problem. The parameters are the same as in the previous example. To 
compare to a single-porosity simulation, we consider a smaller 5 m x 5 m domain with three 



discrete fractures, as indicated in Figure [14| for CF elements with a large width of 10 cm. The 
fracture aperture is 4 mm. Gas, with the same composition as in the previous example, is injected 
from the bottom at 10% PV/yr and production is at constant pressure from the top. 



Figure [15] shows the three phase compositions of C0 2 at 30% PVI. There are small differences, 
but the main flow pattern is captured well by the CF simulation, considering that the CF simulation 
required about 2 mins, while the single-porosity simulation took 74 hrs, a factor of over 2100. The 
average CFL time-step for the CF simulation is about 1.3 day, while the time-steps for the single- 
porosity simulation are about 7 x 10" 4 day, which accounts for most of the difference in total CPU 
time. The remaining improvement in CPU time is achieved by the more efficient pressure update, 
due to the reduced contrast in permeability between the CF and matrix elements. 

4.6. Example 6: Large 3D domain with ten discrete fractures 

In this last example, we consider a large 600 m x 100 m x 50 m three-dimensional domain 
containing ten discrete planar 2D fractures. Each fracture plane is characterized by two diagonally 
opposite corners with coordinates (*i,;yi,Zi) and (jc 2 , yi,Zi), respectively, as provided in Table [6} 
We consider 4 fractures in the y-z direction, and 3 fractures in both the x-y and x-z orientations. 
The domain is discretized by 48 x 29 x 14 elements with 1 m width of the CF elements. The 
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mesh and locations of the fractures, injection and production wells are illustrated in Figure [16} We 
consider about six orders of magnitude range in spatial scales, and four orders of magnitude in 
permeability: the fractures have an aperture of 1 mm and permeability of 100 d, while the matrix 
blocks have a permeability of 4 md and porosity of 44%. 

The domain is initially saturated with the fluid characterized in Table [7] at a temperature of 
400 K and a pressure of 300 bar at the bottom. In the simulation, one pore volume of C0 2 is 
injected at a rate of 10% PV/yr from the corner atOmxOmxOm and production is from 



the opposite corner at 600 m x 100 m x 50 m at constant pressure. Figure [17] shows the overall 
C0 2 composition throughout the domain at 1%, 8%, 25% and 100% PVI. Not surprisingly, we 
find that once the C0 2 front reaches the first fracture, C0 2 quickly flows through the connected 
fractures, resulting in early breakthrough. In this example, we neglected Fickian diffusion, so there 
is no efficient mechanism for cross-flow between the fractures and the neighboring matrix blocks. 
This example demonstrates that we can apply our discrete fracture model to model compositional 
multiphase flow in complex three-dimensional discretely fractured domains. 

5. Summary and conclusions 

First, we briefly reiterate the advantages of our higher-order finite element modeling of three- 
phase flow, and the discrete fracture model. The main motivation for the use of the combined 
MHFE and DG methods is flow in heterogeneous and fractured porous media. From MHFE, we 
obtain an accurate pressure field across interfaces of difference permeabilities (fractures), because 
the pressures are continuous across element edges as well as inside the elements through the use of 
the Raviart-Thomas basis vector fields. At the same time, a continuous velocity field is provided 
at every point, and to the same order of convergence. This is a marked improvement over lowest- 
order FD models that only update average pressures in each element, and compute velocities as 
a post-process. Phase properties, on the other hand, are intrinsically discontinuous across per- 
meability jumps. The DG method is therefore particularly suitable for the mass transport update, 
rather than approximating discontinuous properties with continuous methods, even at higher order. 
Specifically, at edges between a fracture and a matrix block, we have a single continuous pressure 
and corresponding flux, but a different composition in the fracture from the matrix. Within each 
element, we use a higher-order approximation to the mass-transport update, which reduces the 
numerical dispersion as compared to lower-order methods. The higher-order convergence means 
that accurate results can be obtained on coarser grids and at lower CPU cost. 

We adopt the CFE approximation to model fractures, which has considerable advantages over 
commonly used single- and dual-porosity models. Compared to single-porosity models, the CPU 
cost is reduced by orders of magnitude by combining fractures with a small amount of matrix fluid 
into larger computational elements, which relaxes the CFL condition. The CPU time is further im- 
proved because the CF elements have a lower contrast in (effective) permeability with the matrix 
blocks, which speeds up the matrix inversion in the pressure update. As long as the pore volume 
of the matrix-slice included in the CF elements is small compared to the total pore volume of ma- 
trix blocks, results from the CF approximation are indistinguishable from fine mesh simulations. 
Compared to dual-porosity models, the CF approach has the advantage that there is no need for 
transfer functions, which may not have a solid physical footing. Interactions between fracture and 
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matrix elements are computed as in single-porosity simulations. Challenging examples that benefit 
from this approach, and often require discretized matrix blocks, include gravitational and viscous 
fingering in matrix blocks, gravitational re-infiltration of oil from fractures to neighboring matrix 
blocks, and Fickian diffusion. The sugar-cube configuration, required by dual-porosity models, 
may also be overly restrictive, for example in modeling discrete fractures. 
In this work, we have advanced these methods in a number of areas: 

• We have presented the first discrete fracture simulator for fully compositional EOS -based 
three-phase flow, including both three hydrocarbon phases with transfer of all species be- 
tween the three phase, and two hydrocarbon phases and an aqueous phase with C0 2 solubil- 
ity in water. 

• We have adopted a CFL condition for compositional multi-phase flow in a fractional flow 
formulation. The corresponding time-step selection greatly improves the numerical stability 
of the method compared to earlier work, and in some cases increases the CPU efficiency. 

• Fickian diffusion is modeled with a self-consistent model based on a full matrix of com- 
position dependent diffusion coefficients, derived from irreversible thermodynamics. The 
open-space diffusion coefficients are reduced by the formation porosity (and tortuosity) to 
account for the reduced area available for diffusion. Diffusion from fracture to matrix ele- 
ments is improved, with respect to single-porosity models, by mixing the fracture fluid with 
a small amount of matrix fluid, which allows fracture-matrix diffusion within the oil and 
gas phases. Further improvements can be made by enforcing thermodynamic equilibrium 
across the fracture-matrix interface, which we consider in a future work. The generalization 



of this work to account for capillarity is presented in Moortgat and Firoozabadi (2013) 



• The detailed description of the MHFE implementation of the CF discrete fracture model 
includes terms omitted in earlier work, in particular related to the gravitational flux in the 
fractures. This improves the results when gravity is an important driving force, such as dur- 
ing depletion. The formulation is also significantly simplified by introducing a weighted 
effective mobility for the CF elements that takes into account both the fracture and the ma- 
trix contributions to the flux. In terms of this mobility, the implementation on both unfrac- 
tured and fractured domains is nearly identical. Our main challenge in modeling fractured 
domains, particularly in 3D, has been to develop an appropriate mesh generator that also ini- 
tializes the simulations by working out all the fracture orientations and the intersection areas 
of fractures with the faces of CF elements. Once these geometric factors are initialized, the 
modeling of fractured reservoirs is surprisingly straightforward with our proposed model. 
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A. Coefficients in MHFE expansion of convective fluxes 

The Pk,e,e' coefficients in Eq. ( fT2| ) for matrix elements are given (as in |Moortgat and Firooz^ 



abadil ( [20T0l ))by: 



— if.-*- 



n-i 



(31) 



Eq. pTj ) expresses that for each edge/face E,0k,e,e is the effective mobility times the length/surface 
(2D/3D) of edge/face E divided by the distance from the mid-point of E to the center of K (denoted 
by L E ±). As an example: on a 2D rectangular mesh, with scalar absolute permeability K m , and 
after performing mass-lumping onto the diagonal, Eq. pi) reduces to: 

fi K ,E,E^ =2k t ^ K ^l h + l fl v Y (32) 

where we have defined the diagonal matrixes I h = diag[l, 1,0,0], I v = diag[0, 0, 1, 1] and l x 
and l y are the lengths of a rectangular element in the x and y directions, respectively (with edges 
numbered right = 1, left = 2, top = 3, bottom = 4). 

Similarly, for horizontal and vertical fractures (with s the area of the fracture intersection with 
E) 9 we have in 2D: 

/? M ^(hor) = 2ek ftK I h /l X9 /? M ^(vert) = 2sk LK I v /l y . (33) 
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For a cross-flow element we sum the contributions from the matrix (fik,E,E r ) and fracture ifif, e , e f ) 
portions. For example, a 2D element K containing a horizontal fracture has 

l y — s s 

Pkx\ = A,u + PfXi = k f , m — — -I/* + — Ittjlh, and /?^ 2 ,2 = Pkxi- ( 34 ) 

In terms of the total weighted effective mobility across edge E, as defined in Eq. ( [T5] ), we can 
succinctly write for the mass-lumped Pk,e,e' 



@K,E,E' - /^K,E,E^E,E f ~ Pk,E,E f + Pf,e,e f ~ k| 



eff A £ 



1 L 



(35) 



-£± 



in both 2D and 3D, and for both matrix elements and fracture-containing CF elements. In terms 
of k eff , the remaining coefficients in Eq. ( [T3] ) reduce to the same form as in Eq. ( [T2| ): 



@K,E - Ok,E + Qf,e - /.0K,E,E', 

E' 

7k,e = 7k,E + 7f,e = - ^Pa,K^f K (g • n^A^. 



(36a) 
(36b) 



B. Coefficients in MHFE expansion of the pressure equation 

The coefficients in Eq. ( [T8| ) are defined in terms of the weighted effective fractional flow func- 
tions f^ K as: 

WUK - /. c a,K^ia,Kfa,K^ an( ^ S ',^ = /_^ C a,KXia,Kfa,K^*a,K' (37) 



and in Eq. ( [T9| ), we have 



$i,K,E - I S^ 



• n KtE . 



E'edK 



JK 



EGdK 



( n c 

VK,E7K,E + 2_j Vi,KSi,K,L 
i=l 



(38) 



w,x£ is as in Eq. (37), but may be evaluated at the edges from the DG results. We define 

n c 
V K ,E = ^ Vi,Km,K,E, (39) 

and write the coefficients in Eq. ( [20] ) as: 

a K = ^Vk,e0k,e, (40a) 

EeK 

Pk,e = J] Vk,ePk,e,e', (40b) 



(40c) 
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C. Matrices in global MHFE velocity and pressure systems 



The matrices R and M, and vector I in Eq. ( [23] ) collect the coefficients defined in Appendix |A| 
for each element and edge: 



ReR n k ,n e ? R KE = e K , E , 

Mel" , M E9E ,= Yj Pw> 



leR NE , I 



K:E,E'edK 

E = X 7k,e- 

KiEedK 



(41a) 
(41b) 

(41c) 



The matrices in the global system for the pressures, Eq. (J25J), are defined similarly from the coef- 
ficients in Appendix [Bj 



DeR M , D KK = j- + a K , 

At 

ReK M , R K , E =h,E 

€PK(t yi) 



G6H * . 0,-^^+^+^ F.-ZZ«iE« 



£ Of 



(42a) 
(42b) 

(42c) 
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Table 1: Initial composition (mole fraction) z?, acentric factor oj, critical temperature T c , critical pressure p c , molar 
weight M w , critical volume V c and volume translation s for the fluid characterization in Example 1. 



Species 




to 


T C (K) 


/> c (bar) 


M w (g/mole) 


V c (cm 3 /g) 


s 


H 2 


0.00 


0.344 


647 


221 


18 


2.14 


0.000 


C0 2 


0.00 


0.239 


304 


74 


44 


2.14 


0.020 


Ci-Nj 


0.45 


0.011 


190 


46 


16 


6.14 


-0.154 


c 2 -c 3 


0.12 


0.118 


328 


47 


35 


4.73 


-0.095 


C4-C6 


0.07 


0.234 


458 


34 


70 


4.32 


-0.047 


C6-C9 


0.08 


0.370 


566 


26 


108 


4.24 


0.038 


C10-C15 


0.12 


0.595 


651 


19 


166 


4.31 


0.115 


co 16+ 


0.16 


1.427 


824 


10 


386 


3.75 


0.277 



Table 2: Density and viscosity for water, oil and CO2 at p = 300 bar and T = 350 K in Example 1. 



Density (g/cm 3 ) 
Viscosity (cp) 



HoO 



Oil 



CO? 



0.985 
0.36 



0.713 
0.53 



0.754 
0.03 



Table 3: Relative permeability parameters for all numerical examples. 





c0 
^ rw 


s° 

^ row 


c0 
°rog 


c0 

°rg 


k° 

^rw 


k° 

^TOW 


k° 


k° 


«w 


^OW 


n og 


ft g 


Ex. 1 


0.0 


0.5 


0.1 


0.0 


0.3 


1.0 


0.6 


1.0 


3.0 


3.0 


3.5 


2.4 


Exs. 2-3 


0.3 


0.5 


0.1 


0.0 


0.3 


1.0 


0.4 


0.6 


3.0 


2.0 


2.0 


2.0 


Exs. 4-5 


0.05 


0.2 


0.2 


0.05 


0.65 


0.5 


0.5 


0.65 


3.0 


3.0 


3.0 


3.0 


Ex.6 














1 


1 


1 


1 


1 


1 


1 


1 
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Table 4: Initial composition (mole fraction) z?, acentric factor co, critical temperature T c , critical pressure p c , molar 
weight M W9 critical volume V c and volume translation s for fluid characterization in Example 2. 



Species 




CO 


T C (K) 


p c (bar) 


M w (g/mole) 


V c (cm 3 /g) 


s 


H 2 


0.00 


0.344 


647 


221 


18 




2.14 


0.000 


co 2 


0.00 


0.239 


304 


74 


44 




2.14 


0.100 


Ci-Nj 


0.57 


0.012 


189 


46 


16 




6.09 


-0.157 


c 2 -c 3 


0.16 


0.120 


330 


46 


35 




4.73 


-0.094 


C4-C6 


0.08 


0.233 


455 


35 


69 




4.32 


-0.048 


C6-C10 


0.09 


0.428 


584 


24 


120 




4.25 


0.055 


co 11+ 


0.11 


1.062 


751 


13 


293 




4.10 


0.130 



Table 5: Density and viscosity for water, oil and CO2 at p = 300 and 350 bar and T = 350 K in Example 2 and at 
p = 350 bar and T = 400 K in Example 3. 





H 2 


Oil 


co 2 


Gas 


Density (g/cm 3 ), 
Viscosity (cp), 


p = 350 bar, T -- 
p = 350 bar, T -. 


= 350K 
= 350K 


0.987 
0.37 


0.586 
0.20 


0.841 
0.04 


— 


Density (g/cm 3 ), 
Viscosity (cp), 


p = 300 bar, T -. 
p = 300 bar, T -. 


= 350K 
= 350K 


0.985 
0.37 


0.602 
0.23 


0.782 
0.03 


0.260 
0.03 


Density (g/cm 3 ), 
Viscosity (cp), 


p = 300 bar, T -. 
p = 300 bar, T -. 


= 400K 
= 400K 


0.953 
0.22 


0.567 
0.16 


0.560 
0.03 


0.237 
0.03 



Table 6: Fracture characterization in Example 6. 



Index 


xi (m) 


y\ (m) 


Z\ (m) 


x 2 (m) 


?2(m) 


Z2(m) 


1 


100 


5 


5 


100 


80 


45 


2 


240 


25 


25 


240 


75 


40 


3 


410 





10 


410 


80 


50 


4 


560 


10 


5 


560 


90 


45 


5 


300 


25 


5 


500 


25 


35 


6 


100 


50 


25 


250 


50 


45 


7 


10 


75 


25 


150 


75 


45 


8 


450 





10 


600 


50 


10 


9 


10 


10 


20 


300 


30 


20 


10 


250 


30 


40 


500 


80 


40 
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Table 7: Initial composition (mole fraction) z?, acentric factor oj, critical temperature T c , critical pressure p c , molar 
weight M w , and volume translation s for fluid characterization in Example 6 (we note that the volume shift for C35+ 
may be lower than for a real reservoir oil). 



Species 




CO 


T C (K) 


Pc(bar) 


M w (g/mole) 


s 


co 2 


0.0083 


0.239 


304. 


74. 


44. 


0.020 


Q-N 2 


0.4427 


0.011 


190. 


46. 


16. 


-0.154 


c 2 -c 3 


0.1177 


0.118 


328. 


47. 


35. 


-0.095 


C4-C6 


0.0741 


0.234 


458. 


34. 


70. 


-0.047 


C7-C9 


0.0821 


0.370 


566. 


26. 


108. 


0.038 


C10-C15 


0.1158 


0.595 


651. 


19. 


166. 


0.115 


C16-C22 


0.0551 


0.870 


717. 


16. 


247. 


0.169 


C23 _ C34 


0.0465 


1.060 


797. 


14. 


336. 


0.223 


C35+ 


0.0577 


1.100 


958. 


10. 


558. 


0.010 



• = Pk cell averaged pressure; o = tpx,E edge averaged pressure 




fine mesh elements 



cross-flow equilibrium element 



Figure 1 : Illustration of cross-flow equilibrium approximation. 
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Figure 2: Example 1 — Computational mesh: 2 m x 10 m domain, 5 cm wide CF fracture elements. Grid 1 has 
1 1 x 57 elements, and grid 2 is 21 x 105. Injection and production wells are indicated by ellipses near the bottom and 
top, respectively. 
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Figure 3: Example 1 — Oil saturation at 100% PV water injection for single-porosity simulation on Grid 1, and CF 
simulations on Grids 1 and 2. 
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Figure 4: Example 1 — Oil saturation at 100% PV water injection and 250% PV CO2 injection for single-porosity 
simulation on Grid 1, and CF simulations on Grids 1 and 2. 
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Figure 5: Example 1 — Oil recovery. 
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Figure 6: Example 2 — Computational grids for 50 m x 50 m domain, 5 m x 5 m matrix blocks, 20 cm wide CF 
fracture elements. Injection and production wells are indicated by circles. 
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Figure 7: Example 2 — Gas (left) and oil (right) saturation after depleting 5% PV in one year on Grid 3. 
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Figure 8: Example 2 — Water saturation at 5% (left) and 50% (right) PV of water-flooding on Grid 3. 
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Figure 9: Example 2 — Overall, gas, oil and water molar composition of C0 2 at 100% PVI of C0 2 on Grid 3. 
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Figure 10: Example 2 — Oil recovery. 
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Figure 1 1 : Example 3 — CO2 molar composition in aqueous phase (left) and aqueous phase density in kg/m 3 (right) 
at 30% PV of C0 2 injection, with K m = 100 md on Grid 2. 
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Figure 12: Example 4 — Overall, gas, oil and intermediate-phase molar composition of C0 2 at 100% PVI on Grid 2. 
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1 3 : Example 4 — Oil recovery. 
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Figure 14: Example 5 — Computational mesh: 5 m x 5 m domain, 10 cm wide CF fracture elements, 24 x 24 grid 
elements. Injection and production wells are indicated by circles in the bottom-left and top-right corners, respectively. 
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Figure 15: Example 5 — Gas, oil and intermediate-phase molar composition of CO2 at 30% PVI. 
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Injector 




Producer 



Figure 16: Example 6 — Computational mesh: 600 m x 100 m x 50 m domain, 1 m wide CF fracture elements, 
48 x 29 x 14 grid elements. An iso- surface plot of the 50 d permeability level gives an indication of the locations of 
the 10 discrete fractures. 
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Figure 17: Example 6 — Iso-surface plots of the 25 mol% level for the overall CO2 concentration at 1%, 8% and 
25% PVI, and contour plot for the overall CO2 concentration at 100% PVI. 
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